|
|
† Corresponding author. E-mail:
Project supported by the National Natural Science Foundation of China (Grant No. 61372172).
The additional sparse prior of images has been the subject of much research in problems of sparse-view computed tomography (CT) reconstruction. A method employing the image gradient sparsity is often used to reduce the sampling rate and is shown to remove the unwanted artifacts while preserve sharp edges, but may cause blocky or patchy artifacts. To eliminate this drawback, we propose a novel sparsity exploitation-based model for CT image reconstruction. In the presented model, the sparse representation and sparsity exploitation of both gradient and nonlocal gradient are investigated. The new model is shown to offer the potential for better results by introducing a similarity prior information of the image structure. Then, an effective alternating direction minimization algorithm is developed to optimize the objective function with a robust convergence result. Qualitative and quantitative evaluations have been carried out both on the simulation and real data in terms of accuracy and resolution properties. The results indicate that the proposed method can be applied for achieving better image-quality potential with the theoretically expected detailed feature preservation.
X-ray computed tomography (CT) has had a revolutionary effect in fields ranging from biology to medicine. Given that excessive x-ray radiation exposure may cause genetic disease, the reduction of radiation dose during x-ray CT inspection has been a focus of recent investigations.[1–3] Consequently, strategies that lower the required number of projections have received great attention.[4] Owing to insufficient angular sampling, the reconstruction results usually suffers aliasing artifacts, and conventional reconstruction methods such as the well-known filtered back projection (FBP), the algebraic reconstruction technique (ART),[5] simultaneous ART (SART),[6] and expectation maximization (EM) methods[7] do not reliably converge on the correct solution or afford satisfactory image quality. Reconstructing a target image from few-view data has become a research goal in medical imaging.
Inspired by the compressive sensing (CS) theory, methods involving total variation (TV) minimization have been developed and widely studied for reconstruction from few-view data.[4,8–13] The TV norm is the ℓ1 norm of the image gradient magnitude, and minimizing the ℓ1 norm is known to exploit sparsity in its argument.[4] However, owing to the piecewise constant assumption, TV regularization reconstruction methods may lead to a staircasing effect, which may suffer from over-smoothness in fine structures and produce undesirable images.[14,15] To achieve more accurate results, improved edge-preserving methods that involve the introduction of a penalty weight on the original TV norm with more local information have been proposed.[16–20]
The conventional TV norm uses only the sparsity information of the local continuity of an image, an approach that may encounter limitations in some richly detailed or piecewise polynomial situations. In recent years, a nonlocal (NL) gradient characterizing the structural similarity has been studied and has shown good performance for detail preservation in many image-processing applications.[21–23] A method using the ℓ1 norm of the NL gradient, termed nonlocal total variation (NLTV) as a regularization term, has been proposed to improve image quality.[24–26]
Though NLTV shows potential advantages, unfortunately, it often suffers problems in the application of CT image reconstruction. The calculation of an NL gradient is based on the statistics of similarity in a reference image which reflect the similarity interaction of the reconstructed object. However, this reference is difficult to obtain in sparse-view CT image reconstruction as the observed CT projection data do not directly reflect the similarity information of the object. Thus, in CT image iterative reconstruction, the available reference used to calculate similarity weight is an intermediate result that struggles to approximate the actual result. Compared to the true image, it always contains many fakes and blurs. When calculating the NL operator from these references, it may introduce much false information and lead to an incorrect result. Reference [25] used the FBP image as a reference to build the weight. This method requires all angular projections to guarantee the accuracy of similarity estimation. Reference [26] updated the similarity weight in each iteration to generate a more reliable result. However, because the calculation of the similarity weight needs lots of time and memory, this method will cost a great amount of computation.
This paper intends to develop a robust and efficient reconstruction algorithm to make the most of the sparsity brought from the similarity prior. A hybrid optimization model utilizing both TV and NLTV regularizations is proposed. To solve the optimization problem, we then present an efficient iterative algorithm based on the alternating direction method (ADM) optimization approach.[27–29] The new method shows a good performance in sparsity exploitation of both local continuity and nonlocal spatial self-similarity. In particular, it only needs to take several times of calculation of the similarity weight. Thus, the proposed method could preserve a good balance between performance and computation. Finally, we performed qualitative and quantitative evaluations to verify the feasibility of artifact reduction and detail recovery for sparse-view CT image reconstruction. For simplicity, the present method is termed “HTV-ADM” in this paper.
The motivation of NL methods is to exploit the structural similarity property. In this case, the NL gradient is calculated with pixels belonging to the whole or semilocal region of an image, instead of only the nearest neighboring pixels as used in the conventional gradient model. As described in Ref. [23], the NL gradient of an image u: Ω → R is defined as
Let Uω(x) denote the gray-level vector of a patch Px centered at point x and let Ωω(x) denote a neighborhood around point x. The similarity weight ω(x,y) can be written as follows:[22]
In this work, for a general discussion of sparsity exploring and performance evaluation of the NL gradient, we used a typical weight expression that may show basic similarity affinity for most applications. The following similarity function is examined:
Figure
 | Fig. 1. Calculation of the image NL gradient. For a target pixel x in red, all pixels in its search window (orange square) are matched with x by calculating similarity weights. For example, the similarity between pixel y and x is measured by the weight ω(x,y), and the weight is calculated and depends on an ℓ2 distance between the patches (green squares) centered on the pixels in the comparison. |
In most medical and other tomographic imaging applications, the object image generally has the property of spatial self-similarity. There are often many similar structures in different regions of the image, no matter whether the structure is smooth or complex. An NL gradient penalizes the differences between these similar regions, and thus can exploit the considerable information of the structural similarity prior. This approach will obviously lead to a sparse representation of an image u in the basis ∇NL.
For CT image reconstruction, the similarity weight function 
CT projection and desired image can be mathematically assumed to satisfy the following discrete linear equations:
Inspired by the studies of standard TV in image reconstruction, we specify an HTV model that implements the following optimization problem:
The introduction of the TV and NLTV terms in this optimization process differentiates infinitely many solutions to Eq. (
By introducing vectors
The solution to minimizing 𝓛A can be approximated efficiently by the alternating direction method: at each iteration find an approximate minimizer with respect to the variables
The first-order necessary conditions for optimality are
The exact minimizer of Eq. (
(iv) Finally, the Lagrange multipliers are updated as follows:
The calculation of similarity weight of NL operator usually takes a huge amount of computation and is time-consuming. Based on our experience, the weight does not need to be calculated and updated in every iteration. In our implementation, the similarity weight was set to be updated at regular intervals with a number of sub-iterations to approximate the true results, and it only needs to update two or three times of the weight for most image reconstruction procedure. As the available information in iterative procedure does not yield a good estimation for the similarity weight ω at the start, we consider a weight updating strategy: in the beginning of the iterative procedure, z and rz are not updated and set to zero, and the algorithm is executed the same as the TV-ADM algorithm; when the changes of images between two adjacent iterations become small, then the current image is set as the initial reference, and the weight ω is fixed and updated during reconstruction iterations.
In summary, the implementation of the present HTV-ADM method for x-ray CT image reconstruction can be described as follows:
In the above implementation scheme, ΩNL denotes the set of the sub-iteration number for updating similarity weight and WNL (
For any reconstruction design, multiple parameter settings are likely to be involved, which can have a significant impact on reconstruction results. This section discusses the parameter settings of the NL gradient for guiding adequate adaptation of the specific CT imaging task.
As expressed in Section 2.1, it is easy to see that the similarity weight ω is the kernel of the NL gradient and is affected by four important parameters: i) the size of the patch, ii) the size of the search window, iii) the number of neighbors in Ωω(x), and iv) the similarity scaling parameter h. The details of these parameter settings are described below.
The patch size controls the areas considered to match with each other. In general, a greater patch size means more pixels in each patch. As it is hard to match every pixel between two large patches, a large size will introduce limited information on similarity redundancy. In contrast, a small patch size may be preferred. The complexity of the object image is also an important factor. If the image contains too many intricate details, it will be impossible to obtain sufficient information to discern subtle structures without a small patch size. Based on our experience, a value between 3 and 13 is adequate for most images.
In theory, when the patch size and neighbor number are fixed, a search window of increased size can explore more regions and may update pixels considered as neighbors in the set Ωω(x), as they have higher similarity with the target pixel x. If the size of the search window is too small, the performance of the NL gradient will be similar to that of an anisotropic gradient.
However, increasing the size of the search window also results in a heavy computational demand. To preserve a balance between performance and computation, a realistic approach is to make an appropriate selection according to image content. Based on our experience, a value which is twice as large as the patch size is appropriate.
Neighbor number m plays a key role in specifying the level of similarity. Employing more neighbors could directly increase similarity information, but at the same time, may increase some unsuitable connections that introduce an unfavorable prior. Thus, it appears to seek a balance point beyond which similarity information used increases little, while the influence of unsuitable information continues to increase. Based on our experience, a value between 7 and 13 seems to be an appropriate choice for most cases. Another feasible choice is to use a moderate threshold to limit the weight value.
The scale parameter h acts as a soft threshold controlling the decay of the exponential function and determines the norm values considered similar. In general, h often corresponds to the noise level, and is fixed as the standard deviation of the noise.[24] Thus, it is necessary to build a noise estimation model for choosing h.
As a single scale parameter h may not be optimal for all the patches in the image,[31] in this paper we use an adaptive local estimation inspired by local scaling model[32]
The advantage of localizing the parameter h is that it can control the difference penalty according to the region’s content by studying the local statistics of the neighborhoods.
To demonstrate and validate our new method for the CT image reconstruction, two types of data (computer-simulated digital CS-phantom projection data, and experimental rabbit head projection data) were used. All of the experiments were performed under Matlab 2012a running on an HP-Z820 workstation with Intel Xeon E5-2650 dual-core CPU 2.60 GHz. To bring forth the advantages of the presented HTV-ADM algorithm, the TV-ADM[28] algorithm was also performed for comparison.
In the simulation study, a modified CS-phantom[33] designed for compressed sensing reconstruction was used to simulate the few-view projection data. For reference, the used CS-phantom was shown in Fig.
 | Fig. 2. (a) A modified CS-phantom shown in the grayscale window [0,1]. (b) A gradient magnitude image of the phantom shown in the grayscale window [0,0.05]. The phantom array comprises 52780 nonzero pixels, and there are 22779 nonzero values in its gradient image. |
 | Table 1. Parameters in simulation studies. . |
In the first group of simulation, the number of iterations for all methods was 2500. For TV-ADM and HTV-ADM algorithms, the parameters of ADM framework were the same in the experiments: μ and λ1 were set to 1024 and 32, respectively. For the HTV-ADM algorithm, α1, α2, and λ2 were set to 1, 1, and 32, respectively. The NL gradient was updated at iterations 500 and 1000, and the patch size and searching window were set to 9×9 and 31×31, respectively.
The images reconstructed from ideal data by FBP, TV-ADM, and HTV-ADM methods are shown in Fig.
 | Fig. 3. Image reconstruction of the CS-phantom from noise-free projection dataset. Results obtained from (a) FBP, (b) TV-ADM, and (c) HTV-ADM methods. Display window is [0, 1]. |
 | Fig. 4. Reconstructed ROIs of the CS-phantom from noise-free projection dataset. Results obtained from (a) original image, results of the (b) TV-ADM, and (c) HTV-ADM methods. Display window is [0.3, 0.6]. |
To evaluate the reconstruction accuracy quantitatively, the root mean squared error (RMSE) is used as a measure of the reconstruction deviation between the reconstructed and ideal images. The RMSE is defined as
The corresponding RMSE of the reconstructions is calculated and plotted in Fig.
 | Fig. 5. Root mean squared errors as a function of iterations for different methods in noise-free case. |
| Table 2. RMSE and running time of different methods in the noise-free case. . |
In the simulation of a noisy case, admissible reconstructions need more projection views than the noiseless case because of inconsistencies in data acquisition. Thus, simulations with noisy data were conducted using 120 views. For the TV-ADM and HTV-ADM algorithms, μ and λ1 were set to 128 and 32, respectively. For the HTV-ADM algorithm, α1, α2, and λ2 were set to 1, 1, and 32, respectively. As the noise would easily be introduced into the reconstructed image with the increase in iteration number, the iteration number for TV-ADM was selected as 200. The NL gradient was updated at iterations 200, and the total number of iterations for HTV-ADM was 500. The patch size and searching window were set to 11×11 and 31×31, respectively.
The images reconstructed by the FBP, TV-ADM, and HTV-ADM methods are shown in Fig.
 | Fig. 6. Image reconstruction of the CS-phantom from noisy projection dataset. Results obtained from (a) FBP, (b) TV-ADM (iteration number: 200), and (c) HTV-ADM (iteration number: 500) methods. Display window is [0, 1]. |
 | Fig. 7. Reconstructed ROIs of the CS-phantom from noisy projection dataset. Results obtained from (a) original image, (b) TV-ADM (iteration number: 200), and (c) HTV-ADM (iteration number: 500) methods. Display window is [0.3, 0.6]. |
 | Fig. 8. Root mean squared errors as a function of iterations for different methods in the noisy case. |
| Table 3. RMSE and running time of different methods in the noisy case. . |
To further demonstrate the potential capability of the proposed method, we performed an experimental rabbit head study from real data for clinical use. In this study, CT projection data were acquired using a CT scanner primarily comprising an x-ray source (FXE-225, YXLON, Germany) and a flat-panel detector (Varian 4030E, USA). The central slice of the projection data was extracted for 2D investigation and was modeled with 1600 bins on a 1D detector for 2D image reconstruction. The associated imaging parameters of the CT scanner were as follows: (i) 360 projection views were acquired evenly for a 360° rotation on a circular orbit, (ii) the distance from the detector to the x-ray source was 1434.73 mm, (iii) the distance from the rotation center to the source was 502.808 mm, and (iv) the space of each detector bin was 0.254 mm. All of the reconstructed images comprised a 720×720 square of pixels. The size of each pixel was 0.089 mm × 0.089 mm. We evenly extracted a 72-, 120-, and 180-view projection from the projection data for illustration purposes.
In the experiment, the algorithm parameters of all ADM frameworks were set to μ = 128 and λ1 = 64. For the HTV-ADM algorithm, α1, α2, and λ2 were set to 1, 1, and 64, respectively. Similar to Section 3.1.2, the numbers of iterations for the TV-ADM and HTV-ADM were 200 and 600, respectively. The NL gradient was updated at iterations 200, and the patch size and searching window were set to 11×11 and 31×31, respectively.
The reconstructed images of the rabbit head from the real CT scan are presented in Fig.
 | Fig. 9. Images reconstructed using the (a) FBP, (b) TV-ADM, and (c) HTV-ADM methods from 360- (1st row), 180- (2nd row), 120- (3rd row), and 72-view (4th row) projections, respectively. Display window is [0.0047, 0.08] mm−1. |
 | Fig. 10. Zoomed-in views of images reconstructed using the (a) FBP, (b) TV-ADM, and (c) HTV-ADM methods from 360- (1st row), 180- (2nd row), 120- (3rd row), and 72-view (4th row) projections, respectively. Display window is [0.0047, 0.08] mm−1. |
Sparsity in the image gradient is exploited widely for reducing the x-ray CT sampling rate in the projection view angle, but the sparsity in the NL gradient is embedded and often ignored in state-of-the-art CT image reconstruction. To extract this information, we investigate the sparse representation of the to-be-estimated image by an NL gradient, and then improve the optimization-based reconstruction by exploring the local continuity and nonlocal similarity prior. Experiments have shown a clear improvement over standard TV minimization algorithms owing to the introduction of the self-similarity information of spatial interactions into the image reconstruction process.
As the sparsity of the NL gradient is related to the exploitation of image similarity redundancy, the setting of the similarity description plays a vital role and has considerable value to be explored and discussed. For different parameter settings, the reconstruction is likely to yield different levels of image quality. In this study, we focus on the sparsity-guided metrics of parameter setting determinations to use most of both similarity and sparsity potential for a CT imaging model. Although we cannot provide a “best” selection strategy for this HTV minimization procedure, the results show a similar performance in accordance with our analysis. The finding that the suggested method in this study employing HTV minimization allows accurate image recovery with sparse projection data suggests a clinically useful potential for radiation dose reduction.
The framework and metrics are only considered for the typical l2 distance-based similarity weight. Other affinity measures between regions and pixels may provide a more reasonable performance. The complexity of the algorithm is also an important issue. As the computational complexity of weight computing is linear, it will increase rapidly with increasing image size. We hope to extend the theoretical foundations and also to investigate effective affinity expressions for gaining much improvement at very low computational cost. Addressing this question is one of our future research focuses.
| 1 | |
| 2 | |
| 3 | |
| 4 | |
| 5 | |
| 6 | |
| 7 | |
| 8 | |
| 9 | |
| 10 | |
| 11 | |
| 12 | |
| 13 | |
| 14 | |
| 15 | |
| 16 | |
| 17 | |
| 18 | |
| 19 | |
| 20 | |
| 21 | |
| 22 | |
| 23 | |
| 24 | |
| 25 | |
| 26 | |
| 27 | |
| 28 | |
| 29 | |
| 30 | |
| 31 | |
| 32 | |
| 33 | |
| 34 |
